Comput.  Methods  Appl.  Mech.  Engrg.  217-220  (2012)  46-57 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Comput.  Methods  Appl.  Mech.  Engrg. 

journal  homepage:  www.elsevier.com/locate/cma 


A  3D  interface-enriched  generalized  finite  element  method  for  weakly 
discontinuous  problems  with  complex  internal  geometries 

Soheil  Soghrati3,  Philippe  H.  Geubelle5* 

a  Department  of  Civil  and  Environmental  Engineering,  University  of  Illinois  at  Urbana-Champaign,  205  North  Mathews  Avenue,  Urbana,  II  61801,  USA 
b  Department  of  Aerospace  Engineering,  University  of  Illinois  at  Urbana-Champaign,  104  South  Wright  Street,  Urbana,  IL  61801,  USA 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  24  July  2011 

Received  in  revised  form  14  November  201 1 

Accepted  27  December  201 1 

Available  online  4  January  2012 


Keywords: 

GFEM/XFEM 
Weak  discontinuity 
Convection-diffusion  equation 
Heat  transfer  problem 
Convergence  study 
Heterogeneous  materials 


An  interface-enriched  generalized  finite  element  method  (GFEM)  is  introduced  for  3D  problems  with 
discontinuous  gradient  fields.  The  proposed  method  differs  from  conventional  GFEM  by  assigning  the 
generalized  degrees  of  freedom  to  the  interface  nodes,  i.e.,  nodes  generated  along  the  interface  when  cre¬ 
ating  integration  subdomains,  instead  of  the  nodes  of  the  original  mesh.  A  linear  combination  of  the 
Lagrangian  shape  functions  in  these  integration  subelements  are  then  used  as  the  enrichment  functions 
to  capture  the  discontinuity  in  the  gradient  field.  This  approach  provides  a  great  flexibility  for  evaluating 
the  enrichment  functions,  including  for  cases  where  elements  are  intersected  with  multiple  interfaces. 
We  show  that  the  method  achieves  optimal  rate  of  convergence  with  meshes  which  do  not  conform  to 
the  phase  interfaces  at  a  computational  cost  similar  to  or  lower  than  that  of  conventional  GFEM.  The 
potential  of  the  method  is  demonstrated  by  solving  several  heat  transfer  problems  with  discontinuous 
gradient  field  encountered  in  particulate  and  fiber-reinforced  composites  and  in  actively-cooled  micro- 
vascular  materials. 
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1.  Introduction 

The  generalized/extended  finite  element  methods  (GFEM / 
XFEM),  first  introduced  in  [1-4],  is  now  widely  used  by  the  compu¬ 
tational  mechanics  community  to  solve  problems  with  strong  or 
weak  discontinuities,  i.e.,  problems  presenting  discontinuities  in 
the  solution  or  gradient  fields,  respectively.  By  eliminating  the 
need  for  meshes  that  conform  to  the  internal  geometry  of  the  prob¬ 
lem  and/or  propagating  discontinuities  independently  of  the  finite 
element  discretization,  these  methods  provide  a  level  of  flexibility 
absent  in  standard  finite  element  methods  (FEM).  This  unique 
capability  is  achieved  by  incorporating  a  priori  information  on 
the  solution  field  in  the  numerical  approximation  with  the  aid  of 
enrichment  functions.  For  problems  with  weak  or  strong  disconti¬ 
nuities,  appropriate  enrichment  functions  are  selected  to  capture 
the  jump  in  the  gradient  or  solution  fields,  respectively,  allowing 
for  meshes  that  are  independent  of  the  morphology  of  the  prob¬ 
lem.  Application  areas  of  the  GFEM/XFEM  include  fracture 
mechanics  [5-10],  multiscale  modeling  [11],  contact  problems 
[12,13],  solidification  [14],  modeling  of  dislocations  [15,16],  and 
phase  interfaces  [17-19].  A  comprehensive  review  of  the  GFEM/ 
XFEM  method,  the  associated  enrichment  functions,  and  their 
implementation  issues  is  provided  in  [20]. 
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In  the  current  paper,  we  present  an  interface -enriched  general¬ 
ized  finite  element  method  (IGFEM)  for  3D  problems  with  C°- 
continuous  fields  such  as  those  encountered  in  thermal  and 
structural  problems  involving  composites  [21-23],  polycrystalline 
materials  [24,25],  fluid/solid  interaction  [26],  and  conjugate  prob¬ 
lems  with  multiple  phases  [27].  Unlike  in  conventional  GFEM 
where  the  generalized  degrees  of  freedom  (dofs)  are  assigned  to 
the  nodes  of  the  original  mesh,  the  IGFEM  assigns  the  generalized 
dofs  to  interface  nodes  created  by  the  intersection  of  the  phase 
interface  with  the  non-conforming  mesh.  Armero  and  co-workers 
[28]  have  proposed  the  introduction  of  element-based  interface 
enrichments  for  strongly  discontinuous  problems  where  the  gener¬ 
alized  dofs  are  associated  with  the  elements  themselves.  However, 
the  approach  adopted  in  the  present  work  relies  on  the  enrichment 
of  the  interface  nodes  and  thereby  enforcing  continuity  of  the 
enriched  solution  between  adjacent  elements. 

A  unique  feature  of  the  IGFEM  is  the  evaluation  of  the 
enrichment  functions,  which  are  obtained  through  a  linear  combi¬ 
nation  of  the  standard  Lagrangian  shape  functions  of  the  integra¬ 
tion  elements.  Hence,  regardless  of  the  orientation  and  number 
of  interfaces  intersecting  with  the  edges  of  an  element,  the  IGFEM 
provides  a  simple  approach  to  construct  the  enrichment  functions 
and  can  be  considered  as  an  h-hierarchical  method  [29].  In  this 
work,  we  present  a  strategy  for  creating  the  integration  elements 
and  the  corresponding  enrichment  functions  for  tetrahedral  ele¬ 
ments.  The  formulation  of  the  IGFEM  for  2D  problems  has  been 
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previously  presented  in  [30],  together  with  a  convergence  study 
and  a  detailed  comparison  of  the  this  method  with  conventional 
GFEM.  In  the  current  work,  besides  extending  the  IGFEM  formula¬ 
tion  to  3D  problems,  we  put  more  emphasis  on  how  the  IGFEM 
handles  geometrical  complexities  involved  in  these  problems  and 
study  the  application  of  this  method  for  solving  several  engineer¬ 
ing  problems  with  intricate  internal  geometries. 

The  outline  of  this  paper  is  as  follows:  In  the  next  section,  we 
present  the  strong  and  weak  forms  of  the  convection-diffusion 
equation  describing  the  thermal  response  of  actively  cooled 
microvascular  composites  that  motivated  this  work.  It  should  be 
noted,  however,  that  although  the  thermal  convection-diffusion 
equation  is  the  focus  of  this  work,  the  IGFEM  can  readily  extended 
for  other  problems,  including  those  characterized  by  strong  discon¬ 
tinuities  as  briefly  discussed  later  in  the  manuscript.  In  the  remain¬ 
der  of  that  section,  we  review  some  of  the  key  features  of 
conventional  GFEM  approximations  before  introducing  the  3D 
IGFEM  formulation  based  on  a  discretization  with  tetrahedral  ele¬ 
ments.  A  short  discussion  is  also  provided  in  that  section  to  com¬ 
pare  the  computational  cost  of  the  3D  IGFEM  with  that  of 
conventional  GFEM.  We  then  analyze  in  Section  3  the  performance 
of  the  IGFEM  by  performing  a  detailed  convergence  study  with 
non-conforming  meshes  and  comparing  the  results  with  those  of 
the  standard  FEM  based  on  matching  meshes.  Finally,  we  present 
in  Section  4  the  IGFEM  solution  to  several  heat  transfer  problems 
with  weak  discontinuities  to  illustrate  the  capabilities  of  the  meth¬ 
od  to  handle  thermal  problems  with  a  variety  of  intricate  internal 
geometries. 

2.  Convection-diffusion  equation  and  GFEM/IGFEM 
formulations 

Consider  the  open  domain  Q  =  Qs  u  Qf  c  IR3  and  its  closure  Q. 
The  domain  is  partitioned  into  two  separate  regions  Qs  and  Qf  cor¬ 
responding  to  solid  and  fluid  phases,  respectively.  A  velocity  field 
v  :  Q  — >  DR3  is  then  defined  for  the  convective  term,  with  v  =  0  in 
Qs  where  only  conduction  contributes  to  the  heat  transfer.  The 
boundary  r  =  Q  -  Q  with  the  outward  unit  normal  vector  n  has 
been  divided  into  three  mutually  exclusive  partitions  ru,  rq ,  and 
rh.  Here,  ru  corresponds  to  Dirichlet  boundary  conditions  with 
the  prescribed  temperature  u  :  ru  — >  IR,  rq  corresponds  to  Neu¬ 
mann  boundary  conditions  with  the  heat  flux  q  :  rq  — ►  IR,  and  rh 
represents  the  region  with  Robin  (convective)  boundary  conditions 
with  the  heat  transfer  coefficient  h  :  rh  — ►  R,  and  ambient  temper¬ 
ature  Uoo  :  rh  — >  IR.  Assuming  the  thermal  conductivity 
k  :  Q  — ►  U3  x  R3,  density  p  :  Q  — ►  IR,  and  specific  heat  cp  :  Q  — ►  IR 
as  the  material  properties  of  the  domain  and  given  the  heat  source 
/  :  Q  — ►  R,  the  strong  form  of  the  convection-diffusion  equations  is 
then  expressed  as:  find  u  :  Q  — ►  IR  such  that: 

-  V  •  (fcVu)  +  pcpv  •  Vu  =/  on  O, 
u  =  ii  on  ru, 

(1) 

KVu  n  =  q  on 

icVu  •  n  =  h^Uoo  -  u)  on  rh. 

The  weak  form  of  (1)  is  expressed  as  follows:  assuming  the  solution 
field  to  belong  to  the  function  space  U  such  that 
U  =  {u  :  u\ru  =  u}  c  H1^)  and  given  the  function  space 
W  =  {w  :  w\Fu  =  0}  c  H^(Q)  as  the  space  of  the  weight  functions, 
find  u  eU  such  that  for  all  w  e  W: 

/  Vw  •  (kVu)  dQ  +  pcpwv  ■  Vu  dQ  +  /  hwu  drh 

Jq  jQf  Jrh 

=  /  wf  dQ  +  /  wqdrq,+  [  hwu^  drh.  (2) 

Jq  Jrq  Jrh 


2.1.  Conventional  GFEM  formulation 

The  GFEM  formulation  is  based  on  the  Galerkin  method  where 
Uh  c  U  and  VhcV  are  selected  such  that  Uh  =  Vh.  The  domain  Q 
is  then  discretized  into  M  finite  elements,  Q  9*  Qh  =  u^O,-,  which 
do  not  need  to  conform  to  the  phase  interfaces  in  Q.  The  standard 
Lagrangian  shape  functions,  {N,-  (x)  :  x  — >  IR}"=1 ,  are  used  to 
approximate  the  solution  field  (and  thus  the  weight  functions) 
in  elements  not  cut  by  any  interface,  where  n  is  the  number  of 
nodes  in  each  element.  To  retrieve  the  missing  information  in 
the  solution  field  in  elements  located  over  the  interface,  a  set  of 
m  enrichment  functions  at  each  node  of  the  element, 
{(pij{x)  :  x  — >  is  employed  to  evaluate  the  gradient  disconti¬ 

nuity.  The  GFEM  approximation  of  the  solution  field  for  (2)  is 
then  expressed  as: 

n  n  m 

uh(x )  =  ^N,'(x)uf  +  ^N,(X)  22  (Pij(x)Uij.  (3) 

i= 1  i=  1  j= 1 

The  first  term  in  (3)  represents  the  standard  FEM  part  of  the 
approximation  although  u,  may  not  represent  the  nodal  values 
of  the  solution  if  the  enrichment  functions  cpy  do  not  vanish  at 
the  nodes.  The  second  term,  i.e.,  the  contribution  of  the  enrich¬ 
ment  functions  to  the  solution  field,  uses  the  concept  of  the 
partition  of  unity  =  1)  to  localize  the  effect  of  enrich¬ 

ment  functions  in  each  element  and  provide  a  uniform  enrichment 
along  the  interface.  It  must  be  noted  that,  for  certain  enrichment 
functions,  the  enriched  nodes  might  be  extended  to  elements  adja¬ 
cent  to  those  intersected  by  the  interface,  known  as  blending  ele¬ 
ments,  to  maintain  optimal  rate  of  convergence  [31].  This 
correction  is  required  since  the  partition  of  unity  is  unable  to 
reproduce  fully  the  enrichment  functions  in  the  blending  elements 
and  thus  corrective  enrichment  functions  are  employed  to  create  a 
smooth  transition  between  standard  FEM  elements  and  those  cut 
by  the  interface. 

Although  the  implementation  of  the  GFEM  is  very  similar  to 
that  of  the  standard  FEM,  two  GFEM-related  implementation  is¬ 
sues  need  to  be  addressed.  The  first  one  is  the  quadrature  in  the  en¬ 
riched  elements  due  to  weak  or  strong  discontinuities  in  the 
enrichment  functions.  The  most  common  approach  to  achieve  an 
accurate  quadrature  in  these  elements  consists  in  dividing  them 
into  subdomains,  referred  to  as  integration  sub-elements,  with 
edges  that  conform  to  the  phase  interfaces  passing  through  the  ele¬ 
ment  [5,6].  The  aspect  ratio  of  these  integration  elements  does  not 
affect  the  accuracy  of  the  solution  as  they  only  serve  as  subdo¬ 
mains  for  choosing  appropriate  integration  points.  The  second  is¬ 
sue  appears  at  enriched  nodes  located  over  the  boundaries  with 
Dirichlet  boundary  conditions,  i.e.,  prescribed  values  of  the  tem¬ 
perature  in  the  current  work.  Since  the  solution  field  at  each  node 
is  obtained  from  the  contribution  of  both  the  standard  and  general¬ 
ized  dofs,  and  since  the  enrichment  functions  do  not  necessarily 
disappear  at  the  boundary  nodes,  the  prescribed  field  values  can¬ 
not  be  directly  assigned  at  such  nodes.  In  this  case,  other  ap¬ 
proaches  such  as  Lagrange  multipliers  or  penalty  method  must 
be  employed  to  enforce  the  prescribed  values  of  the  solution  field 
at  enriched  nodes  [32,33]. 

2.2.  IGFEM  formulation 

The  major  difference  between  the  IGFEM  and  conventional 
GFEM  consists  in  relocating  the  generalized  dofs  from  nodes  of 
the  original  mesh  to  those  of  integration  sub-elements  located  over 
the  interface,  referred  to  as  interface  nodes.  We  then  assign  an 
enrichment  function  to  each  interface  node,  remove  the  partition 
of  unity  from  the  evaluation  of  the  enrichment  functions,  and  ex¬ 
press  the  IGFEM  approximation  as: 
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uh(x)  =  jrNj(x)Gi  +  (4) 

i=  1  i=l 

As  before,  the  first  term  in  (4)  corresponds  to  the  standard  FEM  por¬ 
tion  of  the  approximation.  Since  the  generalized  dofs  are  attached 
to  the  interface  nodes,  the  standard  dofs  in  this  first  term  always 
represent  the  field  values  at  the  nodes  of  the  original  mesh.  The  sec¬ 
ond  term  in  (4)  is  the  contribution  of  the  enrichment  functions  in 
approximating  the  solution  field  along  the  phase  interface  where 
ujnt  are  the  generalized  dofs.  The  enrichment  functions  in  IGFEM, 
i.e.,  (i/q(x)  :  x  — >  R}"=1,  must  capture  the  gradient  jump  in  the 
solution  field  inside  an  element  while  vanishing  at  the  nodes  and 
edges  of  that  element  which  are  not  cut  by  the  interface.  The  latter 
condition  for  the  enrichment  functions  is  necessary  to  avoid  discon¬ 
tinuities  alongs  the  edges  of  the  enriched  elements  and  adjacent 
elements,  i.e.,  to  ensure  the  conformity  of  the  approximate  field 
along  the  boundary  of  the  elements,  and  thus  avoid  any  corrections 
for  blending  elements.  As  explained  below,  the  parameter  s  is  a 
scaling  factor  used  to  avoid  sharp  gradients  in  the  enrichment 
functions. 

In  this  work,  we  describe  the  evaluation  of  the  enrichment  func¬ 
tions  for  tetrahedral  elements,  although  the  approach  can  be  read¬ 
ily  extended  to  other  types  of  elements.  Fig.  1  illustrates  four 
different  scenarios  by  which  a  four-node  tetrahedral  element 
may  be  cut  by  an  interface  and  the  associated  division  of  the  ele¬ 
ment  into  integration  subdomains.  The  approach  adopted  here 
consists  in  creating  the  minimum  number  of  integration  elements 
for  accurate  quadrature  over  the  parent  (original)  element  to 
facilitate  the  evaluation  of  enrichment  functions.  A  tetrahedral  ele¬ 


Fig.  1.  Evaluation  of  enrichment  functions  for  the  3D  IGFEM:  four  scenarios  of 
creating  the  integration  subdomains  composed  of  tetrahedral  and  wedge  elements, 
and  the  corresponding  numbering  of  integration  elements  used  to  evaluate  the 
enrichment  functions. 


ment  is  then  divided  into  a  combination  of  tetrahedral  and 
pentahedral  (wedge)  elements  based  on  the  orientation  of  the 
interface  with  respect  to  the  element  edges.  To  capture  the  gradi¬ 
ent  discontinuity  along  the  interface,  the  enrichment  functions 
corresponding  to  Fig.  l(a)-(d)  are  expressed  as: 

(a)  ifo  =  N™  +  N?\  <A2  =  AfW  +  Nf,  =  N $■>  +  N®, 

<1*4  =  +  Nf] , 

(b)  ^1=N',)  +N'2),  ,A2=N<,)+N<2),  i/,3=N<1)+N<2),  (5) 

(c)  ^  =  N[v  +  N<2) ,  ^  =  N™  +  Nf  +  Nf , 

(d)  f,  =N<1)+N<2), 


respectively.  Please  note  that  the  creation  of  wedge  element  in  the 
IGFEM  is  not  required  and  that  a  parent  element  can  be  discretized 
using  tetrahedral  elements  only,  for  which  the  enrichment  func¬ 
tions  are  constructed  differently  based  on  the  number  of  integration 
sub-elements. 

A  possible  issue  with  the  aforementioned  approach  for  evaluat¬ 
ing  the  enrichment  functions  may  occur  when  integration  ele¬ 
ments  with  very  high  aspect  ratios  are  created,  i.e.,  when  the 
interface  intersecting  the  tetrahedral  elements  happens  to  be  very 
close  to  one  of  the  nodes  of  the  original  mesh.  Since  the  Lagrangian 
shape  functions  of  these  sub-elements  contribute  to  the  evaluation 
of  the  enrichment  functions  in  the  parent  element,  the  high  aspect 
ratio  leads  to  high  gradient  values  and  thereby  to  an  ill- 
conditioned  stiffness  matrix.  To  address  this  issue,  we  employ  a 
scaling  factor  s  as  introduced  in  (4)  to  avoid  excessively  large 
gradients  in  the  enrichment  functions  [19].  To  evaluate  this  param¬ 
eter,  we  adopt  the  approach  suggested  in  [30]  and  define  s  as: 


/2  min(xi,x2)\2 

V  *,+x2  J  ’ 


(6) 


where  Xi  and  Xi  are  shown  in  Fig.  2. 

Eq.  (6)  represent  a  parabolic  function  defined  over  the  edge  of 
the  parent  element  holding  an  interface  node,  with  zero  value  at 
the  defining  nodes  of  this  edge  and  unity  in  the  middle  of  it.  As 
the  interface  approaches  one  of  the  nodes  of  the  element,  the 
scaled  enrichment  function  goes  to  zero  and  vanishes  when  the 
interface  coincides  with  an  FEM  node.  It  must  be  emphasized  that 
the  scaling  factor  s  has  a  constant  value  in  an  element  cut  by  the 
interface  and  the  quadratic  function  introduced  in  (6)  has  only 
been  implemented  to  evaluate  this  value. 

A  key  advantage  of  the  IGFEM  is  its  ability  to  handle  problems 
with  complex  geometries  or  microstructures  by  creating 
enrichment  functions  for  special  cases.  Ref.  [30]  describes  the  eval¬ 
uation  of  enrichment  functions  for  the  case  where  two  or  more 
interfaces  intersect  inside  an  element  for  2D  problems.  This  is 
achieved  by  creating  integration  elements  that  conform  to  the 
geometry  of  the  intersecting  interfaces  inside  the  parent  element 
and  using  the  linear  interpolation  of  the  corresponding  Lagrangian 


Fig.  2.  Geometrical  parameters  *i  and  Xjused  in  the  scaling  factor  entering  the 
enrichment  function  associated  with  interface  node  1  (Eq.  (6)). 


S.  Soghrati,  P.H.  Geubelle /  Comput.  Methods  Appl.  Mech.  Engrg.  217-220  (2012)  46-57 


49 


(8)  can  then  simulate  the  discontinuity  in  the  solution  field  along 
the  interface. 

3.  Convergence  study 

Two  example  problems  are  presented  in  this  section  to  study 
the  convergence  rates  of  the  3D  IGFEM.  In  these  examples,  we 
use  as  reference  solution  either  a  closed-form  solution  or  the  solu¬ 
tion  obtained  with  the  standard  FEM  with  a  highly  refined  con¬ 
forming  mesh.  The  IGFEM  results  obtained  with  non-conforming 
meshes  are  then  compared  with  those  provided  by  the  standard 


Fig.  3.  Creation  of  integration  elements  for  the  evaluation  of  the  enrichment 
functions  in  the  IGFEM  for  an  element  cut  by  two  interfaces. 


shape  functions  as  the  enrichment  functions.  The  same  approach 
can  be  extended  to  3D  problems  where  a  combination  of  tetrahe¬ 
dral  and  wedge  integration  elements  that  conform  to  the  phase 
interfaces  are  used  to  discretize  the  parent  element  and  evaluate 
the  associated  enrichment  functions.  Similarly,  the  IGFEM  provides 
the  framework  for  the  evaluation  of  enrichment  functions  for  more 
complex  cases,  where  most  conventional  GFEM  formulations  tend 
to  fail,  such  as  cases  where  elements  are  intersected  by  multiple 
interfaces.  Fig.  3  shows  one  of  the  possible  scenarios  of  a  four-node 
tetrahedral  element  cut  by  two  interfaces  and  the  corresponding 
integration  elements  created  to  evaluate  the  enrichment  functions, 
which  are  obtained  as: 

h  =N<2)+Nf,  f2=N^+Nf\  ^3=N<2)+fV<3), 

l^4  =  N<1)+N‘2),  iA5  =  +  Nf\  ^6  =  N(31)+N?).  (7) 

The  enrichment  function  for  other  orientations  of  the  element  with 
respect  to  the  phase  interfaces  can  be  evaluated  similarly. 

Many  of  the  comments  made  in  [30]  to  compare  the  2D 
implementation  of  the  GFEM  and  IGFEM  apply  to  the  3D  case 
discussed  here.  For  instance,  enforcing  the  Dirichlet  boundary  con¬ 
ditions  in  the  IGFEM  is  similar  to  that  of  the  standard  FEM  at  the 
nodes  of  the  elements  intersected  by  the  interface  because  the  gen¬ 
eralized  dofs  are  moved  to  the  interface  nodes.  Also,  the  computa¬ 
tional  cost  of  the  3D  IGFEM  with  tetrahedral  elements  compares 
favorably  with  that  of  the  conventional  GFEM,  especially  for  cases 
where  blending  elements  are  introduced.  It  can  be  shown  that  the 
number  of  generalized  dofs  introduced  at  the  interface  nodes  in  the 
IGFEM  is  equal  to  that  of  conventional  GFEM  with  no  correction  in 
the  blinding  elements.  Moreover,  considerably  more  generalized 
dofs  are  added  in  conventional  GFEM  if  enriching  all  the  nodes  of 
the  blending  elements  is  necessary  for  achieving  the  optimal  rate 
of  convergence.  In  other  words,  the  cost  of  the  numerical  solution 
via  IGFEM  for  3D  problems  is  equal  to  that  of  the  best  case  of  con¬ 
ventional  GFEM. 

Finally,  although  we  only  address  in  the  current  work  problems 
with  weak  discontinuities,  constructing  the  enrichment  functions 
for  problems  with  strong  discontinuities  within  the  IGFEM 
framework  is  straightforward.  The  enrichment  functions  needed 
for  this  class  of  problems  can  simply  be  obtained  by  switching 
the  sign  used  in  the  linear  combination  of  the  Lagrangian  functions 
in  (5).  For  instance,  the  enrichment  functions  associated  with  the 
integration  elements  shown  in  Fig.  1  for  problems  with  strong  dis¬ 
continuities  can  be  written  as: 

•A,  =N'1)-N‘2),  V'2=n<,)-n<2),  *3  =  N<1)-N‘2).  (8) 


(C) 


Fig.  4.  (a)  Problem  description,  and  (b)  GFEM  solution  for  the  first  example 
The  opposite  signs  Of  the  Lagrangian  shape  functions  in  each  inte-  problem,  (c)  Temperature  profile  along  the  z-axis  obtained  with  the  IGFEM  solver 

gration  element  used  for  evaluating  the  enrichment  functions  in  for  both  the  convergence  study  and  the  patch  test. 
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(a)  (b) 

Fig.  5.  (a)  Structured  conforming  mesh  used  for  the  standard  FEM  solution,  and  (b)  unstructured  non-conforming  mesh  used  for  the  IGFEM  solution  of  the  first  example 
problem  investigated  in  the  convergence  study. 


FEM  with  conforming  meshes  with  the  same  level  of  refinement  to 
verify  the  optimal  rate  of  convergence.  We  provide  the  conver¬ 
gence  results  in  the  L2-  and  H1- norm  of  the  error,  expressed  as: 


3.1.  Example  1:  Constant  gradient  jump  over  a  flat  interface 

The  dimensions  and  boundary  conditions  of  the  first  example 
problem  are  presented  in  Fig.  4.  The  cubic  domain  is  composed 
of  two  regions  along  the  z-axis  with  conductivity  values  of  1  and 
8  W/m  K  for  the  lower  and  upper  sections,  respectively.  The 
boundary  conditions  for  this  problem  are  prescribed  temperatures 
5  and  10  °C  along  the  top  and  bottom  surfaces,  respectively,  while 
the  other  surfaces  are  insulated.  With  the  origin  of  the  Cartesian 
coordinates  system  located  at  the  lower  front  corner  of  the  cubic 
domain,  a  distributed  heat  source /(x,y,z)  =  5000(z2  -  z  +  1)  is  ap¬ 
plied,  yielding  the  following  exact  temperature  field  in  the  domain: 

-33750Z4  +  833. 3z3  -  2500z2  +  70.9z+ 10  if  z <  0.05 
-4218. 7Z4  +  104.2z3  -312.5z2  +8.86Z  +  7.56  ifz^0.05. 

(10) 

Fig.  4  illustrates  the  IGFEM  solution  field  for  this  problem  with  a 
nonconforming  mesh  built  over  al6xl6xl6  grid.  As  expected, 
the  thermal  conductivity  mismatch  generates  a  jump  in  the 
gradient  field  at  z  =  5  cm.  This  gradient  discontinuity  is  accurately 
captured  by  the  IGFEM  as  shown  in  Fig.  4,  which  depicts  the  tem¬ 
perature  profile  along  the  z-axis. 

Before  the  convergence  study,  we  present  a  modified  patch  test 
by  removing  the  heat  source  applied  over  the  domain  and  showing 
how  the  IGFEM  can  capture  the  gradient  discontinuity  along  the 
interface.  The  exact  solution  for  this  problem  is  a  linear  variation 
of  the  temperature  along  the  z-direction,  with  a  slope  discontinuity 
at  the  interface,  as  illustrated  with  the  solid  green  line  in  Fig.  4.  The 
same  figure  also  provides  the  values  obtained  at  the  nodes  of  the 
non-conforming  mesh  (circles)  and  the  temperature  value  ob¬ 
tained  right  at  the  interface  (square  symbol),  as  provided  by  the 
IGFEM  enrichments.  The  perfect  match  between  the  analytical 
and  numerical  results  in  this  simple  example  shows  that  the 
IGFEM  satisfies  the  patch  test,  illustrates  the  temperature  profile 
along  the  z-axis  obtained  from  the  IGFEM  solution  with  the  non¬ 


conforming  mesh  shown  in  Fig.  5(b),  for  this  patch  test.  The  circu¬ 
lar  nodes  depicted  on  this  temperature  profile,  which  is  composed 
of  two  line  segments,  show  the  IGFEM  solution  at  the  FEM  nodes, 
while  the  rectangular  node  depicts  the  interpolation  of  the  temper¬ 
ature  value  at  the  phase  interface  in  a  nonconforming  element  cut 
by  this  interface.  As  shown  there,  the  IGFEM  provides  super  con¬ 
vergent  results  for  this  problem  and  yields  the  exact  solution  value 
along  the  interface  in  a  nonconforming  element. 

The  convergence  study  for  the  original  problem  is  conducted 
using  conforming  structured  meshes  for  the  standard  FEM,  while 
unstructured  non-conforming  meshes  are  used  for  the  IGFEM  solu¬ 
tions  (Fig.  5))  so  that  elements  located  across  the  interface  are 
intersected  along  different  orientations.  Fig.  6  presents  the  conver¬ 
gence  rates  of  the  standard  FEM  and  IGFEM  for  this  problem,  with 
respect  to  both  the  I2-  and  H1- norm  of  the  error.  As  shown  there, 
the  rate  of  convergence  of  the  IGFEM,  without  employing  any  cor¬ 
rections  to  the  blending  elements,  is  equivalent  to  that  of  the  stan¬ 
dard  FEM.  As  expected,  due  to  the  simplicity  of  the  geometry,  the 
accuracy  of  the  standard  FEM  for  the  same  number  of  dofs  is  supe¬ 
rior  to  that  of  the  IGFEM,  although  the  difference  between  the  two 
methods  is  very  small. 

3.2.  Example  2:  Non-constant  gradient  jump  over  a  curved  interface 

In  this  second  example,  we  study  the  convergence  rate  of  the 
IGFEM  for  a  problem  with  a  cylindrically  shaped  material  interface 
as  shown  in  Fig.  7(a).  The  cubic  domain  with  the  length  of  10  cm 
contains  a  cylinder,  extended  along  the  y-axis,  with  the  diameter 
of  7  cm  and  a  conductivity  value  five  times  larger  than  that  of 
the  surrounding  material.  The  prescribed  temperature  along  the 
bottom  surface  of  the  domain  is  ub  =  0  °C  while  a  heat  flux 
q  =  105(x  —  0.05)  +  500  W/m2  is  applied  along  the  top  surface. 
Both  faces  in  the  x  direction  have  conductive  (Robin)  boundary 
conditions  with  a  heat  transfer  coefficient  of  h  =  8  W/m2  K  and 
an  ambient  temperature  u ^  =  20  °C,  while  the  other  two  surfaces 
are  assumed  to  be  insulated.  The  temperature  field  obtained  from 
the  IGFEM  solution  over  the  domain  and  in  half  of  the  matrix  cut  at 
x  =  5  cm  are  shown  in  Fig.  7(b)  and  (c),  respectively. 

The  IGFEM  and  standard  FEM  rates  of  convergence  of  this 
second  example  are  presented  in  Fig.  8.  Since  no  analytical  solu¬ 
tion  is  available  for  this  problem,  we  use  a  standard  FEM  solu¬ 
tion  with  a  very  refined  conforming  mesh  built  on  a 
85  x  85  grid  as  the  reference  solution.  As  in  the  first  example, 
the  standard  FEM  solutions  for  this  problem  are  obtained  using 
conforming  meshes  while  the  unstructured  non-conforming 
meshes  similar  to  those  shown  in  Fig.  5(b)  are  adopted  for  the 
IGFEM  solution. 
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Fig.  6.  Convergence  rates  in  I2-norm  and  F^-norm  of  the  error  with  respect  to  the  mesh  size  (ft)  and  total  number  of  dofs  (N)  for  the  first  example  problem  shown  in  Fig.  4. 
Conforming  and  nonconforming  meshes  similar  to  those  presented  in  Fig.  5  are  used  for  the  standard  FEM  and  IGFEM  solutions,  respectively. 


(b) 


(c) 


Fig.  7.  (a)  Problem  description  and  IGFEM  thermal  solution  over,  (b)  the  entire  domain,  and  (c)  the  matrix  in  half  of  the  domain  cut  at  x  =  5  cm. 
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Fig.  8.  Convergence  rates  in  L2-norm  and  hF-norm  of  the  error  with  respect  to  the  mesh  size  (ft)  and  number  of  dofs  (N)  for  the  second  example  problem  described  in  Fig.  7. 
(e)  presents  the  convergence  rate  values  obtained  from  the  two  most  refined  FEM  meshes  for  each  plot. 


Fig.  9.  First  application  problem,  (a)  Geometry  and  boundary  conditions  for  a  volume  element  of  particulate  composite  with  alumina  inclusions  embedded  in  a  glass  matrix; 
(b)  Distribution  of  embedded  inclusions. 


As  shown  in  Fig.  8,  the  IGFEM  yields  convergence  rates  sim¬ 
ilar  to  that  of  the  standard  FEM.  Moreover,  for  the  same  num¬ 
ber  of  dofs,  the  accuracy  of  the  IGFEM  is  comparable  to  that  of 
the  standard  FEM.  It  must  be  noted  that,  because  of  the  discret¬ 
ization  error  due  to  the  presence  of  the  curved  interface,  nei¬ 


ther  the  standard  FEM  nor  the  IGFEM  yield  the  optimal  rate 
of  convergence.  However,  this  example  shows  that  the  IGFEM 
achieves  a  performance  similar  to  that  of  the  standard  FEM 
while  removing  the  cost  and  complexity  of  creating  conforming 
meshes. 
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Fig.  10.  (a)  Temperature  field  in  the  first  application  problem  presented  in  Fig.  9(a);  (b)  Temperature  profile  along  one  of  the  diagonal  planes;  (c)  Temperature  profile  along 
the  z  =  4  mm  plane,  showing  the  ability  of  the  IGFEM  scheme  to  capture  the  gradient  discontinuity  along  the  particle/matrix  interfaces. 


4.  Applications 

The  four  heat  transfer  problems  described  hereafter  are  inspired 
by  applications  found  in  engineering  and  material  science.  These 
problems  are  all  solved  with  the  IGFEM,  using  structured  meshes 
that  do  not  conform  to  the  complex  microstructure  or  geometry 
of  the  problem. 

4.1.  Application  1:  Particulate  composite 

In  this  first  application,  we  use  the  IGFEM  to  extract  the  effec¬ 
tive  thermal  properties  of  a  particulate  composite,  composed  of  a 
glass  matrix  with  ellipsoid  alumina  inclusions.  As  described  in 
[34],  these  thermal  properties  depend  not  only  on  the  volume  frac¬ 
tion  of  the  fillers,  but  also  on  the  size  and  distribution  of  the  inclu¬ 
sions.  Therefore,  the  evaluation  of  the  effective  properties  of  the 
composite  involves  the  solution  of  several  configurations  of  the 
inclusions  in  the  matrix  [35].  The  complexity  and  cost  of  generat¬ 
ing  conforming  meshes  for  multiple  configurations  render  the  use 
of  the  standard  FEM  very  labor  intensive  and  make  mesh-indepen¬ 
dent  methods  such  as  the  IGFEM  particularly  attractive. 

Fig.  9(a)  illustrates  the  geometry  and  boundary  conditions  of  a 
unit  cell  of  glass/alumina  composite.  A  prescribed  temperature 
ub  =  0°C  and  a  uniform  heat  flux  g  =  100W/m2  compose  the 
boundary  conditions  along  the  bottom  and  top  surfaces,  respec¬ 
tively,  while  other  faces  are  insulated.  The  thermal  conductivity 
values  of  the  glass  matrix  and  alumina  inclusions  are  Kg  =  1.4 
and  Ka  =  31  W/m  K,  respectively.  The  distribution  of  the  alumina 
beads  in  the  glass  matrix  is  presented  in  Fig.  9(b)  where  the 


ellipsoidal  inclusions  have  radii  ranging  from  0.9  to  1.3  mm  and  as¬ 
pect  ratios  between  0.75  and  1.25. 

Fig.  10  illustrates  the  temperature  field  in  this  problem,  ob¬ 
tained  from  the  IGFEM  solution  with  a  non-conforming  structured 
mesh  built  on  a  32  x  32  x  32  grid.  Fig.  10(b)  shows  the  tempera¬ 
ture  field  along  a  diagonal  place,  demonstrating  the  heterogeneous 
nature  of  the  thermal  solution  inside  the  composite.  A  slice  of  the 
temperature  profile  along  the  plane  z  =  4  mm  is  presented  in 
Fig.  10(c),  which  clearly  illustrates  the  ability  of  the  IGFEM  to  cap- 


(a)  (b) 

Fig.  11.  Second  application  problem:  (a)  Geometry  and  applied  boundary  condi¬ 
tions  of  a  unit  cell  of  3D  woven  composite  made  of  epoxy  matrix,  glass  weft  and 
warp  fiber  tows,  and  copper  z-fibers;  (b)  Configuration  of  the  weft  and  warp  tows, 
and  of  the  z-fibers. 
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ture  the  gradient  discontinuity  along  the  perimeter  of  the  alumina  of  the  temperature  field  is  presented  along  the  centerline  of  the  z- 
inclusions.  fiber  at  y  =  9  mm. 


4.2.  Application  2:  3D  woven  composites 

In  the  second  example,  we  use  the  IGFEM  to  evaluate  the  tem¬ 
perature  field  in  a  unit  cell  of  a  3D  woven  composite  inspired  by 
the  material  system  discussed  in  [36].  The  dimensions,  boundary 
conditions,  and  microstructure  of  the  unit  cell  are  shown  in 
Fig.  11.  The  bottom  surface  of  this  50  x  18  x  40  mm  domain  is  sub¬ 
ject  to  a  constant  heat  flux  q  =  400  W/m2,  the  top  surface  has  con¬ 
vective  boundary  conditions  with  h  =  7  W/m2  K  and  u ^  =  20  °C, 
and  other  faces  are  insulated.  The  epoxy  matrix  holds  glass  fibers 
in  the  weft  and  warp  directions,  each  phase  with  thermal  conduc¬ 
tivities  of  Ke  =  0.3  W/m  K  and  Kg  =  0.96  W/m  K,  respectively.  The 
z-fibers  weaving  together  the  weft  and  warp  tows  are  made  of  cop¬ 
per  with  kc  =  401  W/m  K.  A  high-conductivity  material  for  the  z- 
fiber  is  adopted  to  increases  the  thermal  conductivity  of  the  com¬ 
posite  in  the  through-thickness  direction  [37]. 

Fig.  12  illustrates  the  temperature  field  in  the  unit  cell  and  in 
the  fibers.  A  structured  mesh  built  on  a  50  x  24  x  40  grid  is  used 
for  the  evaluation  of  the  IGFEM  solution  for  this  problem. 
Fig.  12(b)  shows  the  role  of  the  copper  z-fiber  in  transferring  the 
heat  through  the  thickness  of  the  unit  cell.  The  figure  illustrates 
the  ability  of  the  IGFEM  to  capture  the  gradient  discontinuity  along 
the  surface  of  the  copper  wire  although  the  conductivity  of  the 
copper  is  more  than  1300  times  larger  than  that  of  the  surrounding 
epoxy.  This  behavior  is  clearly  observed  in  Fig.  12(c)  where  a  slice 


Fig.  13.  Third  application  problem:  (a)  Geometry,  boundary  conditions,  and 
material  properties  of  a  soil  specimen  composed  of  several  layers  where  the 
thermal  conductivity  of  each  layers  is  given  in  W/m  K;  (b)  Temperature  field 
obtained  with  the  IGFEM  solver. 
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Fig.  14.  (a)  Geometry  and  boundary  conditions  for  an  actively-cooled  epoxy  fin;  (b)  Embedded  sinusoidal-shape  microchannels  with  the  amplitude  and  wavelength  of  >4  =  2 
and  A  =  5  mm,  respectively,  and  water  as  the  coolant  with  flow  rates  and  entrance  temperatures  specified  in  the  figure. 


<c) 

40.9 
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(a)  (b) 

Fig.  15.  (a)  Temperature  field  in  the  microvascular  epoxy  fin  shown  in  Fig.  14;  (b)  Temperature  profile  inside  the  domain,  cut  along  the  plane  y  =  1.35  mm,  emphasizing  the 
thermal  impact  of  the  embedded  microchannels. 


4.3.  Application  3:  Laminated  soil 

In  the  third  application,  the  temperature  field  in  a  multi-layered 
soil  sample  with  the  dimensions  shown  in  Fig.  13(a)  is  evaluated 
with  the  IGFEM.  The  thermal  conductivity  of  soil  is  a  function  of 
several  parameters  such  as  the  grain  size,  moisture  content, 
density,  porosity,  and  organic  matter  and  thus  could  be  consider¬ 
ably  different  from  one  layer  to  the  next  [38].  The  temperature 
of  the  specimen  is  fixed  at  ub  =  20  °C  along  the  bottom  surface  of 
the  specimen  and  a  constant  heat  flux  with  intensity  q  = 
1000  W/m  K  is  applied  along  the  top  boundary,  while  other  faces 
are  kept  insulated.  The  IGFEM  solution  obtained  with  a 
25  x  25  x  75  structured  mesh  is  presented  in  Fig.  13(b).  The  tem¬ 
perature  field  is  slightly  warped  in  both  the  x  and  y  directions  to 
emphasize  the  discontinuity  in  the  gradient  field  across  the  inter¬ 
faces  between  different  layers. 

4.4.  Application  4:  Actively-cooled  microvascular  polymeric  fin 

Microvascular  materials  are  being  considered  for  various  engi¬ 
neering  applications  such  as  autonomic  materials  [39,40], 
biotechnology  [41,42],  and  microelectromechanical  systems 
(MEMS)  [43-45].  In  such  materials,  the  diameter  of  the  embedded 
microchannels  varies  from  a  few  microns  to  a  millimeter,  and  the 
embedded  network  of  microchannels  is  designed  based  on  a  vari¬ 
ety  of  application-driven  objective  functions  (thermal  impact,  flow 
efficiency,  mechanical  impact,  . . .)  and  manufacturing  constraints. 
In  this  last  example  problem,  we  employ  the  IGFEM  to  evaluate 
the  temperature  field  in  an  actively-cooled  microvascular  epoxy 


fin  [46-48].  The  specimen  dimensions,  applied  thermal  boundary 
conditions,  and  configuration  of  the  microchannels  are  shown  in 
Fig.  14.  The  10  x  2.7  x  5.2  mm  epoxy  fin  with  k  =  0.19  W/m  K  is 
subject  to  a  constant  heat  flux  q  =  4000  W/m2  along  the  bottom 
surface  and  convective  boundary  condition  with  h  =  7W/mK 
and  Uoo  =  20  °C  along  the  other  surfaces. 

Fig.  14(b)  illustrates  the  configuration  of  the  embedded 
microchannels  in  the  epoxy  matrix  with  diameters  D  =  500  pm. 
A  sinusoidal-shaped  curve  with  the  amplitude  A  =  2  mm  and 
wavelength  2  =  5  mm  describes  the  centerline  of  the  microchan¬ 
nels.  The  coolant  circulating  in  the  channels  is  water  with  kw  = 
0.6  W/m  K,  p  =  1000  kg/m3,  and  cp  =  4185.5  J/kg  K.  The  middle 
channel  carries  a  flow  rate  of  Q.  =  2  ml/min  with  the  entrance  tem¬ 
perature  Te  =  20  °C  while  the  two  other  microchannels  have  a 
slightly  smaller  flow  rate,  Q=  1.5  ml/min,  with  Te  =  15  °C  in  the 
opposite  direction. 

The  evaluation  of  the  convection  term  in  (1 )  associated  with  the 
fluid  flow  requires  the  knowledge  of  the  velocity  field  in  the  micro¬ 
channels.  Due  to  the  very  small  size  of  microchannels  and  the 
small  value  of  the  mass  flow  rate,  laminar  flow  condition  with 
fully-developed  velocity  profile  adequately  describes  the  flow  in 
the  channels.  The  magnitude  of  the  velocity  in  the  flow  direction 
at  a  distance  r  from  the  centerline  of  the  microchannel  is  then  ex¬ 
pressed  as  [49]: 

where  D  is  the  mean  velocity  in  the  channel. 
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The  IGFEM  solution  for  this  problem  is  presented  in  Fig.  15. 
From  Fig.  15a,  which  shows  the  temperature  field  along  three  of 
the  outer  surfaces  of  the  domain,  the  thermal  impact  of  the  coolant 
circulating  in  the  three  microchannels  is  clearly  observable. 
Fig.  15c  presents  details  of  the  temperature  field  inside  the  domain 
along  the  vertical  plane  y  =  1 .35  mm  (at  the  centerline  of  the  mid¬ 
dle  channel),  depicting  the  jump  in  the  gradient  field  along  the  sur¬ 
face  of  the  microchannel  due  to  the  mismatch  in  the  thermal 
conductivity  values. 

5.  Conclusions 

An  interface-enriched  generalized  finite  element  formulation 
has  been  introduced  to  solve  3D  problems  characterized  by  discon¬ 
tinuous  gradient  fields.  In  this  method,  the  generalized  degrees  of 
freedom  introduce  to  enrich  the  numerical  solution  in  non- 
conforming  elements  traversed  by  a  material  interface  are 
associated  with  the  interface  nodes,  and  not  with  the  nodes  of 
the  non-conforming  mesh  as  in  conventional  GFEM.  This  choice 
leads  to  some  implementation  advantages  such  as  the  straightfor¬ 
ward  handling  of  Dirichlet  boundary  conditions.  Moreover,  the 
number  of  the  generalized  dofs  associated  with  the  IGFEM,  which 
determines  the  computational  cost  of  the  method,  is  similar  to  that 
of  conventional  GFEM  in  the  absence  of  blending  elements.  The 
enrichment  functions  for  the  IGFEM  consist  in  a  linear  combination 
of  the  Lagrangian  shape  functions  defined  on  the  integration  ele¬ 
ments.  In  this  work,  the  enrichment  functions  for  four-node  tetra¬ 
hedral  elements  have  been  derived.  It  was  also  pointed  out  how 
the  IGFEM  is  able  to  handle  more  complex  cases  such  as  the  pres¬ 
ence  of  multiple  interfaces  in  a  given  element.  The  flexibility  of  the 
IGFEM  for  evaluating  the  enrichment  functions  makes  this  ap¬ 
proach  highly  suited  for  problems  with  intricate  internal 
geometries. 

A  detailed  convergence  study  of  the  IGFEM  for  3D  heat  transfer 
problems  with  straight  and  curved  interfaces  has  also  been  con¬ 
ducted,  showing  that  the  IGFEM  yields  similar  rate  of  convergence 
and  accuracy  level  with  non-conforming  meshes  as  that  of  the 
standard  FEM  with  conforming  meshes.  Thermal  problems  in  het¬ 
erogeneous  materials  and  structures  ranging  from  composites  to 
microvascular  polymer  have  been  investigated  to  demonstrate 
the  flexibility  of  the  IGFEM  in  capturing  weakly  discontinuous 
solutions  in  a  variety  of  complex  internal  geometries. 
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